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Abstract 

We consider a theoretical model of a four-mode Bose-Hubbard model consisting of two pairs 
of wells coupled via two processes with two different rates. The model is naturally divided into 
two subsystems with strong intra-system coupling and much weaker coupling between the two 
subsystems and has previously been introduced as a model for Josephson heat oscillations by 
Strzys and Anglin [Phys. Rev. A 81, 043616 (2010) ]. We examine the quantum dynamics of this 
model for a range of different initial conditions, in terms of both the number distribution among 
the wells and the quantum statistics. We find that the time evolution is different to that predicted 
by a mean-field model and that this system exhibits a wide range of interesting behaviours. We 
find that the system equilibriates to a maximum entropy state and is thus a useful model for 
quantum thermalisation. As our model may be realised to a good approximation in the laboratory, 
it becomes a candidate for experimental investigation. 

PACS numbers: 03.75.Lm, 03.75.Kk, 67.25.du 
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I. INTRODUCTION 



Bose-Einstein condensates (BECs) of weakly interacting dilute gases have long been recog- 
nised as a valuable tool for the exploration of the dynamics of non-equilibrium many-body 
physics. Experimental investigations of BEC have provided a whole new toolbox for the 
study of quantum mechanics in mesoscopic systems. A recent development is the proposal 
by Strzys and Anglin to use a four-mode Bose-Hubbard model with greatly differing tunnel- 
ing rates as a model for the investigation of mesoscopic thermodynamics, in particular with 
regard to the transport of heat [l^. Their analysis is motivated by the fact that, microscop- 
ically, heat is energy stored in degrees of freedom whose evolution is too quick to perceive 
or control on a macroscopic time scale. The authors performed an analytical analysis of the 



system, based on the treatment given to a two well model by Milburn et al. j2|. In fact, their 
model consists of two of the systems considered by Milburn et al., with a weaker coupling 
between these. 

In this work we will not focus on the analogy Strzys and Anglin make between slow 
Josephson oscillations and second sound jst, but instead will investigate the quantum dy- 
namics and statistics of this system. We find that there is a wealth of complex behaviour, 
very little of which is predicted by mean-field or linearised Bogoliubov type theories. We 
find that, by going beyond linearised analyses and using the fuller, semi-quantum truncated 
Wigner approximation we are able to investigate the relaxation of the system to equi- 
librium. By constructing a measure which behaves very much like von Neumann entropy 
and which is potentially measurable in the laboratory, we show how the four-well model can 
make a contribution to the study of thermalisation in isolated quantum systems jsl . 

II. PHYSICAL MODEL AND HAMILTONIAN 

Extending the standard procedure for two wells we consider a four well potential with 
an independent condensate in each of the four wells at the beginning of our investigations. 
The Hamiltonian for a condensate in an external trapping potential, VJ=a;t(^), may be written 

as 



n= dr 



V fe2 



h ^ ^ „ „ „ „ 

2m 



(1) 
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FIG. 1: Schematic of our four-mode Bose-Hubbard system. The and bi are the bosonic annihi- 
lation operators for each mode, while J and uj represent the coupling rates between the modes. In 
this article, we always set uj = O.IJ and J = 1, which sets the units of time. 

where is the field operator for the condensate, and the non-linear interaction parameter is 
Uq = 2'n'ah/m, where a is the s-wave scattering length describing two-body collisions within 
the condensate, and m is the atomic mass. In the case where the external potential provides 
a four well confinement for the condensate, we may simplify the above Hamiltonian by 
making use of the four-mode approximation. At zero temperature all atoms in the system 
are condensed and if the ground state energies of the condensate in the four single (and 
separate) wells are sufficiently separated from the energies of the condensate in all other 
excited single particle states, transitions to or from the two modes of interest and these 
higher lying states can be neglected. We may then expand the field operator as 

2 
i=l 

where and bi {i = 1,2) are bosonic annihilation operators in each of the wells, and the 
0j are the ground state spatial wave functions of the condensate in wells on the left and 
right side, as shown schematically in Fig. [1] 
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Using this in Eq. ([T]), we find an effective Hamiltonian 

2 

i=l 

+hxa\a\aiai + hxblblbibi^ 
—hJ (^a\a2 + dl^i + b{b2 + blbi 
—hu (^d{bi + blcii + dlb2 + bld2^ , (3) 

where we have neglected the spatial overlap of the different well densities. The single well 
bound state energies, E^^^, are 

E^'' = \jdr (0f/^)*(f) (^V^ + V;.,(f)) 0f/^(f). (4) 

J, the tunnel couphng on each side of the system, is 

^ = ^ / rfr (0[/^)*(r-) (^V^ + K..(r)) ^^/^tf), (5) 

while w, the tunnel coupling between the left and right subsystems, is 

^ = ^jdr (0f )*(r) (^V^ + I4.*(r)) 0f (f). (6) 

The effective non-linear interaction term is 

X = Uoj rff l^f/^ir- (7) 

We set the single well bound state energies equal because we will consider only a symmetric 
potential where we can set = Eji = 0. We note here that, while Strzys and Anglin 
began their dynamical investigations from the ground state and provided a periodic tilt to 
the potentials to excite the dynamics, we leave the potential unperturbed and excite the 
dynamics via differences in the initial populations of the wells. 

We parametrise time by setting J = 1, so that dimensionless time as displayed in the 
results will be in units of Jt. We will investigate the effects of changing x and the initial 
distribution of atoms in the wells. Because we use a quantum phase space method, we may 
also change the quantum statistics of the initial states in each well. In this work we will 
investigate the dynamics arising from initial Fock and coherent states Q]. We will always 
use a total average atom number of A^^ = 10^. 
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III. THEORETICAL METHODS 



For the numerical investigation of many-body interacting quantum systems which are 
too large for master equation methods, the preferred first option is the positive-P represen- 
tation |7|, which allows for an exact mapping from the type of Hamiltonian used here to 
stochastic differential equations. However, in cases where the system is und^amped and has 
a high x^^^ nonlinearity, it tends to become unstable after very short times [8] . As this is the 
case here, we will perform our investigations using stochastic integration in the truncated 
Wigner representation 0, which enables us to capture the majority of the quantum fea- 
tures of the system as long as the Wigner pseudoprobabihty distribution is positive. There 
is no reason to believe this is not the case in any of the investigations we perform here ex- 
cept in the representation of initial Fock states, where we will use an approximation which 
is justified below. The truncated Wigner representation also has the huge operational ad- 
vantage of remaining stable over relatively long integration times. The truncated Wigner 
representation has been shown to be accurate for the investigation of a range of conden- 

nn 

sate dynamics [9Hll| , and we have also found that its predictions are accurate for twin- well 
dynamics, so expect it to be accurate here. 



To find the appropriate equations, we begin by using the operator correspondences [12 1 

«P ^ + 

^ ("* - Wia) (10) 

pat ^ + (11) 

to give a generalised Fokker-Planck equation with third-order derivatives. Although it is 
possible to map this approximately onto stochastic differential equations [IS], the numerical 
integration of these is extremely unstable, so we will instead use what is known as the 
truncated Wigner representation by dropping derivatives of higher than second order in the 
Fokker-Planck equation. This leaves an equation with no diffusion terms for the Wigner 
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pseudoprobability function, 
dW 



dt 



d d 
d d 



0P2 



W. 



(12) 



The equations of motion for the Wigner variables of this system are then found as 

dai 



dt 
da2 

dl3i 
"dt 



—2ix\cii\^ai + iJa2 + 
—21x1(^2^^(^2 + iJai + iujl32 
-2zx|/3iP/3i + iJ(^2 + i^(^i 



^ = -2ix\P2? P2 + iJ Pi + i^(^2, 
dt 



(13) 



where the aj and /3j are the Wigner variables corresponding to dj and 6j, respectively. 
Classical averages of the Wigner variables correspond to symmetrically ordered operator 
expectation values, so that the necessary reordering must be undertaken before we arrive at 
solutions for physical quantities, for which normal ordering is more appropriate. Although 
Eq. (fT3|) might look classical, the Wigner variables themselves are drawn from appropriate 
distributions for the desired initial states, so that the stochasticity comes from the initial 
conditions. The truncated Wigner equations above are solved numerically by taking aver- 
ages over a large number of stochastic trajectories, with initial conditions drawn from the 
distributions given below. 

The dynamical evolution of this system can depend on the initial quantum state as well 
as the initial number distribution. In this case, we will investigate two different initial 
number distributions and two different initial quantum states. Most of our analyses will 
be performed with half the atoms in each of two diagonally opposite wells and the other 
two initially vacant, with equal populations in each well also being used to calculate the 
Josephson frequencies. Using the Wigner representation we may easily simulate different 
initial quantum states [l^. To represent the Wigner distribution for a coherent state, 
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where a\a) = a\a), the initial conditions are chosen from the distribution 

aw = a + ^{i^i + iJ^2), (14) 
where the uj are independent Gaussian normal random variables. We easily see that, as 



required by the symmetric ordering, |avi/P = Na + 1/2. For simulations using initial coherent 
states, we used the open source software package xmds pJi]. Fock states of fixed atom number 
may be simulated using a Gaussian approximation developed by Gardiner et al. [l^ and 
previously used to analyse trapped BEG photoassociation V]\. A Fock state of fixed atom 



number, |iV), can be sampled, to a good approximation as long as N is not too small, by 

aw = {p + qjy)e^'"^ (15) 
where ^ is a random number from the uniform distribution [0, 1) and q = l/4j9, with 

1 / , N 1/2 

p=-[2N + l + 2ViV2 + Nj . (16) 

This approximation has been shown to reproduce well the first two moments for reasonable 
sizes of N 17|, which is all that is required of a Gaussian distribution. Simulations using 
initial Fock states were performed in Matlab which, although not as fast as xmds, does have 
a uniform random number generator. We note here that newer versions of xmds also have 
a uniform random number generator, but this was not available at the time we ran our 
simulations. The initially unoccupied wells have a distribution chosen from Eq. (|T4l) with 
a = 0, which reproduces the vacuum state. 

We note here that we have used these two quantum states not because they are actually 
what we would necessarily expect the physical quantum state of a condensate to be, but 
because they are in common usage and serve to show any dynamical differences which may 
arise from differences in the initial quantum statistics. In fact, other states have been used 
previously, with Olsen and Plimak jol, Q using both squeezed states and states which are 
sheared in the phase-space |18l] to investigate BEG photoassociation. 



IV. QUANTUM AND CLASSICAL DYNAMICS 



When they introduced this system, Strzys and Anglin 
tum djTiamics was rich beyond the scope of their paper 



stated that the full range of quan- 
1|. In this work we reveal a small 



part of this rich dynamics, beginning by exposing the differences from the classical dynam- 
ics as predicted in the coupled Gross-Pitaevskii equation (GPE) approach. In their paper, 
Strzys and Anglin claim that the number of atoms they use (10^) is large enough so that 
the dynamics will stay close to the classical predictions. However, when we use the same 
parameters as in their paper, with J = 1, u = 0.1 and x = 2.5J/Nt, where Nt is the 
total number of atoms, we see that this is not the case. In Fig. |2l we compare the classical 
Gross-Pitaevskii predictions for the populations in wells ai and 02 to those found for initial 
coherent states in the truncated Wigner representation, we see that they are markedly dif- 
ferent after the first few oscillations. We note here that, although we have not used the same 
initial number distribution as Strzys and Anglin, our result shows that having a relatively 
large total number of atoms is not sufficient in itself to mean that the dynamics are a small 
perturbation around the classical predictions. We also note that, due to the remaining sym- 
metry of the system for these parameters and initial conditions, we only need to show the 
populations of one of the two wells on one side to demonstrate the full population dynamics. 

In this work, we will use much lower coUisional nonlinearities (x) than that used by Strzys 
and Anglin, but the dynamics we will investigate will still not be a small perturbation about 
the classical predictions. To illustrate this, as shown in Fig. [3l we take an initial condition 
with half the atoms in each of wells ai and 62 and a nonlinearity of x = O-U/Nt, with the 
two different quantum states, and compare these to the prediction of the classical equations. 
Over the time shown, the prediction for initial coherent states is almost indistinguishable 
from that of the GPE approach, although a collapse in the oscillations is seen at longer 
times. This collapse is normally explained for x^^^ systems as being due to the fact that 
different number components of the coherent state superposition oscillate at different rates, 
so that they eventually become out of phase. Although supersymmetry arguments say that 
atoms cannot actually be in coherent states 19|, we would expect the same dynamics for 
a mixture of number states. The Fock state results, on the other hand, are very different, 
unlike the case of two coupled wells where they only become markedly different for higher 
collisional strengths. 

When we increase the collisional interaction to x = 0.5J/Nt, we see that the mean-field 
dynamics change qualitatively for both initial quantum states, as shown in Fig. |H The 
collapses and partial revivals in the oscillations occur on a time scale inversely proportional 
to the interaction strength and both initial distributions equilibriate to equal numbers in each 
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FIG. 2: (colour online) Populations Nai and Na2 with the atoms initially distributed evenly in 
wells ai and 62. The parameters used are J = 1, w = 0.1 and x = 2.5J/N'r, where Nj- is the total 
number of atoms, with a value of 10^. The dash-dotted lines are the classical predictions, while 
the solid lines are for initial coherent states in the populated wells. The coherent state prediction 
is the result of 6 x 10^ stochastic trajectories. The time axes in this and subsequent time domain 
figures are in dimensionless units. 

of the four wells. For these initial number distributions, and for values of x ^ 2J/A^t, the 
classical prediction is for totally regular oscillations, not displaying any of the macroscopic 



self-trapping predicted in the two-well system {2]. An analysis of the number and quadrature 
statistics shows that the atoms in each well evolve away from being in either coherent or 
Fock states, with more noise in both quadratures and intensity than coherent states, but 
less noise in the quadratures than for Fock states. Another difference is that we did not see 
complete revivals in the oscillations for the four-well system. We will give an explanation for 
their absence in terms of relaxation and the possibility of chaotic behaviour in section IVIIi 
We stress here that this behaviour would not have been observable in the usual linearised 
analyses, which depend on the accuracy of the GPE type approach as a starting point. The 
difference between the classical GPE dynamics and those for the two initial quantum states 
we consider here shows that, for this system, a total number of 10^ atoms is not sufficient 
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FIG. 3: (colour online) Population A^oi with the atoms initially distributed evenly in wells ai and 
62- The parameters used are J = 1, a; = 0.1 and x = 0-1J/Nt, where Nt is the total number of 
atoms, with a value of 10^. The solid line is the classical prediction, with the prediction for initial 
coherent states being identical over this length of time, while the dash-dotted line is for initial 
Fock states in the populated wells. The coherent state prediction is the result of 2 x 10^ stochastic 
trajectories while the Fock state result is the average of 4 x 10^ trajectories. 

for it to be treated classically, but that both the initial quantum state and the quantum 
dynamics must be considered to obtain accurate predictions, even for the mean fields. The 
time over which we have integrated the equations is different for each initial state for three 
reasons. The first of these is that xmds is faster than Matlab, the second is that we need 
more samples to faithfully reproduce Fock states, and the third is that the results with 
initial coherent states generally take longer to demonstrate qualitative differences in their 
dynamics. 

The use of phase-space methods also allows us to calculate the djTiamics of the Schwinger 



pseudospin operators |20| adapted to this four-mode system. Given that there is no diagonal 



tunnelling, that we are examining two linked subsystems, and the symmetry of the initial 
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FIG. 4: (colour online) The longer time population dynamics of wells oi and 02 for the same 
parameters and initial conditions as in Fig. O except that the collisional interaction strength is 
now X = 0.5J/Nt- (a) shows the dynamics for initial coherent states, while [(b)] is for initial Fock 



states. We see that the average number in each well tends to a state where they oscillate about an 
equal distribution, which did not happen for the lower value of x- 

conditions we consider, it is sufficient to define 

Sy = — i ^0.10.2 — a{d2j , 
SP = aSl + a[bi, 

5« = -t ^a,b\ - ajfei) , (17) 

with the obvious changes being made to arrive at Si^^ and Sy^\ As the represent number 
differences, we have not used them here. 

The 5*^ represent the particle occupation number difference between the single party 
energy eigenstates of two adjacent wells, and is clearly related to the coherence between 
these two wells, while the Sy represent the momenta between adjacent wells. We also note 
that we use a different normalisation convention to that used by some authors, who insert 
a factor of l/y/2 in front of the expressions. We will also normalise Si^^ by dividing by 
{a[ai + b\bi) and 5*" by {d\ai + a|a2), so that the range of possible values lies between 
— 1 and 1, with representing no occupation number difference of the eigenstates and ±1 
representing the maximum possible. As can be seen in the Hamiltonian, the Sz are also 
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FIG. 5: (colour online) The longer time dynamics of the momenta between wells ai and 02 for the 



same parameters and initial conditions as in Fig. [3l (a) shows the dynamics for initial coherent 
states, while [(b)] is for initial Fock states. The momenta between bi and 62 are equal in magnitude 
and opposite in sign. 

proportional to the energy of tunnelling between adjacent wells. We will now show that the 
dynamical expectation values of these quantities also depend strongly on the initial quantum 
states. 

Figs. [5] and [6] show the momenta between the two wells within each subsystem, for two 
different collisional nonlinearities and initial quantum states. The differences between the 
two initial states here are markedly qualitative, with obviously different oscillation frequen- 
cies. The longer time behaviour suggests that in all four cases the system tends towards 
an equilibrium situation where the momenta will equal zero, but the envelopes under which 
this behaviour occurs are quite different. We found that the momenta between the two 
subsystems exhibits similar behaviours. While the coherent state collapses and revivals are 
naturally expected, we also see that decaying and not completely periodic behaviour is seen 
for initial Fock states. This behaviour in particular suggests that there are a number of 
frequencies involved here with repeating and almost regular patterns being observable. This 
will be examined further in Section |Vl where we will compare the frequencies found in our 
numerical analysis with those predicted by Strzys and Anglin in their analytic approxima- 
tion, and with those predicted by Bogoliubov theory. 

In Fig. [7] we show Sz between wells ai and 02 for the lower collisional nonlinearity, finding 
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FIG. 6: (colour online) The longer time dynamics of the momenta between wells ai and 02 for 



the same parameters and initial conditions as in Fig. [3l but with x = 0.5J/Nt . (a) shows the 
dynamics for initial coherent states, while [(b)] is for initial Fock states. The momenta between 61 
and 62 are equal in magnitude and opposite in sign. 

once again that the dynamics differ quahtatively for the different initial quantum states. In 
the case of initial coherent states, the oscillations in S"" decay to a finite value, while the Fock 
state case continues to oscillate over the time shown here in a complex, but almost periodic 
manner. We see the same type of behaviour for Si^\ which is between the two halves of the 
system, as shown in Fig. |H]and Fig. |9l What is obvious here is that, when the interaction 
strength is increased, 5*^ and all the other quantities we have calculated settle down faster to 
a type of equilibrium situation where the atoms are evenly distributed among the four wells 
and the momenta and tunnelling energies are not changing very much. In fact, we would 
only expect this system to oscillate indefinitely for a vanishing atomic interaction strength, 
in which case it would be composed of ideal gases and the semiclassical analysis would be 
accurate. The difference in these behaviours is not predictable from the analysis used by 
Strzys and Anglin, but will be seen to be important when we investigate the analogies with 
heat exchange and entropy, below in section IVII 
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FIG. 7: (colour online) The longer time dynamics of the normalised 5" for the same parameters 



and initial conditions as in Fig. [3l (a) shows the dynamics for initial coherent states, while [(b)] is 
for initial Fock states. 
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FIG. 8: (colour online) The longer time dynamics of the normalised si^^ for the same parameters 



and initial conditions as in Fig. [3l (a) shows the dynamics for initial coherent states, while [(b)] is 
for initial Fock states. 



V. JOSEPHSON FREQUENCIES 



Now that we have shown that the classical solutions for the dynamics of this system are 
not always reliable, we will examine our results for evidence of the Josephson oscillations 
predicted in the linear Bogoliubov approximation, as well as the low frequency collective 
mode predicted by Strzys and Anglin. In terms of the units we use, the Josephson frequencies 
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FIG. 9: (colour online) The longer time dynamics of the normalised si^^ for the same parameters 



and initial conditions as in Fig. [3l but with x = 0.5J/Nt- (a) shows the dynamics for initial 
coherent states, while [(b)] is for initial Fock states. 

for elementary excitations above the N-atom ground state are 



cu = 

n' = ^2{J + w)(2J + 2u + Ntx)- 



(18) 



In this case, VL' and VL are at a higher frequency than Co and Stryzs and Anglin also predict 
a beating between these two, with frequency 



lo{2lo + Ntx) 



2u 



4J 



(19) 



when J ^ u. In the cases we examine in this article, J = 1, while a; = 0.1 and NtX has so 
far been either 0.1 or 0.5, so this condition holds. 

We now look for evidence of these frequencies in Fourier transforms of the expectation 
values of the atomic numbers and the S operators, using different initial quantum states 



and number distributions. In Fig. 10(a), we begin with the populations equally distributed 



among the wells in coherent states and take the Fourier transforms of the number in any 
one well. As can clearly be seen, the observed frequencies follow the Bogoliubov predictions 
very closely for this configuration. We note here that the actual number oscillations are 
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FIG. 10: (colour online) The frequencies of the Josephson oscillations found numerically for initial 



coherent states, with the solid lines being the analytical expressions, (a) shows the results for an 
equal initial distribution of atoms in each well, i.e. Naj^{0) = Na^iO) = iV;,i(0) = Nb2i^) = 2,500, 
while [(b)] initially has iVai(O) = A^a2(0) = 2600 and Nb^{0) = Nt^iO) = 2,400 . The tunneling 
interaction strengths are J = 1 and oj = 0.1. 

very small, being driven by the Poissonian number uncertainties in each well. In this case 
we do not see clear evidence in the Fourier components of the beat mode predicted by 
Strzys and Anglin, although we note that it would be only marginally different from u. 
In view of the fact that Strzys and Anglin found their numerical results by adding a time 
dependent potential tilt to each subsystem, we also began with an asymmetric system, with 
2600 atoms in each of the wells on the left hand side, and 2400 in each of the others. As 



shown in Fig. 10(b) , this resulted in one frequency of oscillation only for the atoms, with the 
Fourier transforms being extremely clean, and this frequency closely matching u. We have 
not included frequency analysis of the system beginning in initial Fock states, as we found 
that these needed averaging over too many stochastic trajectories to give clear signals. 

We also examined the frequencies of oscillation in Sy'^^ and Si^\ as these represent tun- 
nelling momentum and interaction between the two subsytems. These are across the weak 
link where Strzys and Anglin expect a process analogous to heat transfer to take place, so 
it is of interest to look for evidence of the collective mode frequencies predicted in their 
paper. When we examine the results for an initial equal number distribution, we find fre- 
quencies closely matching u and in Sy'^\ and two clear frequencies in Si^\ one of which 
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FIG. 11: (colour online) The frequencies of the oscihations found numerically for initial coherent 
states from Sy'^'^ and si^\ with the solid lines being the analytical expressions. The initial atom 



numbers are evenly distributed, (a) shows the results from Sy'^\ while [(b)] is for si^^ . The tunneling 



interaction strengths are J = 1 and u = 0.1. 

matches Cl, while the other is not close to any of the four frequencies we have dealt with so 
far. As far as we can tell, it does not seem to be either of a)± from Eq. (12) of the Strzys 
and Anglin paper, as these should be either comparable to or lower than u, although it is 
difficult to deconstruct their approximations to know what the joson number should be for 
our parameters. These results are shown in Fig. [TTl The frequencies found by beginning 
with the unequal number distribution are shown in Fig. [121 where we again see that not all 
the analytically predicted results are clearly found. We note here that we did not attempt 
to mechanically excite these predicted frequencies, but instead looked for those which may 
arise naturally from population imbalances in the system. 



VI. ANALOGY WITH JOSEPHSON HEAT OSCILLATIONS 

As the principal motivation behind the work of Strzys and Angin was as a contribution 
to the development of a mesoscopic and eventually microscopic model for heat transfer, we 
will now examine the quantities in our system which may be useful in such a model. Heat 
transfer involves two principal processes. The first of these is energy and the second is a 
change in entropy if we begin in a non-equilibrium state, as we do in this article. The 
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FIG. 12: (colour online) The frequencies of the oscihations found numerically for initial coherent 
states from S^y"^ and S^z \ with the solid lines being the analytical expressions. The initial atom 



numbers are Na^{{)) = Na^{Q) = 2600 and iV6i(0) = A^62(0) = 2,400. [(a)] shows the results from 
Sy'\ while [(bj] is for . The tunneling interaction strengths are J = 1 and u) = 0.1. 

expectation values of the various energies involved can be calculated straightforwardly by 
finding the dynamical averages of the quantities in the effective Hamiltonian of Eq. ([3]). As 
expected for a closed system, the total energy is a constant. What is interesting, however, 
is that even as the mean values of the oscillations die down and the atoms approach a state 
where they are evenly distributed among the four wells, the tunnelling energy between sites 



does not disappear, but approaches a constant value. This is shown in Fig. 13(a), where 
we have plotted the different components for the left hand side of the system, where the 
different energies shown are defined as (note that we set ^ = 1 for the graphics so that we 
are not using S.I. units) 

Esite = hx (a\ '^al + al '^afj , 

Ea = E,,te + Etun " ^ (a\hi + h\ai + a\h2 + h\a^ , (20) 

where half the weak link tunnelling energy has been included in the total energy of the left 
hand side. 

We can also calculate a candidate for entropy by following the approach taken by Strzys 
and Anglin in Eq. (13) of their paper. Before we do this, however, we will examine their 
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FIG. 13: (colour online) (a) The different energies in the left subsystem as a function of scaled 
time, for initial Fock states with A^ai(O) = Nt^iO) = 5000 and Na^{0) = Nb^{0) = 0, J = 1, w = 0.1 
and x^T = 0.5. |(b)| The pseudo entropy for the same parameters. The solid line is for initial 
coherent states while the dash-dotted line is for initial Fock states. 

expression for the single-particle reduced density matrix in more detail. Considering only 
one of the subsystems, we may define the even and odd modes excited by a±. = (di ± 02), 
(note that our definition differs by a scale factor of 1/ \/2) which leads to the matrix, 

{a±a±) 



Ra 



{a\di + 62^2) 

1 / (a^a_|_) (a^a_) 



{a\ai + 0^02) I (ala+) (ala_) 



(21) 



which is presented in ref. [1| as a diagonal matrix. However, when we expand their definition, 
we find 

11 {a\ai + a|a2 + S*") {a\ai - 0,20-2 + ^5*^) 

Ra — 1 4^ I I ) (22) 

2(0101 + 0^02) y{a\a^-ala2-iS^) (a{ai + 0^02 + 5^) 
which will not in general be diagonal. Leaving aside the fact that a single-party reduced 
density matrix is only strictly defined in this manner for eigenstates of the number operator 
and we have seen that, even when we begin our system in number eigenstates, it does not 
remain in them, we will follow a similar approach and calculate what we will call a "pseudo 
entropy" for reasons which will become obvious below. 

Rather than using the density matrix defined above, we define an approximate density 
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matrix for the left hand side of the system as 

R = I ^^'^'^ ^^'^'^ I (23) 

where Na = ajo-i + ala2. This allows us to work with the numbers at each of the two sites, 
rather than the numbers in modes which combine both sites. It is then an easy matter to 
calculate a single-particle subsystem pseudo entropy from this matrix, 

e = -Tr(i?Jni?J, (24) 

which will have a maximum value of In 2 ^ 0.6931 when the atoms are equally distributed 
throughout the wells, which is statistically the most probable situation. In Fig. [T3] we 
see that our heuristic entropy does approach this value, although not monotonically. This 
suggests that, while it is formally wrong to think of Eq. ( 123|) as being a density matrix, 
it is only approximately wrong when it comes to calculating the single-particle subsystem 
entropies and could prove to be a useful experimental measure. Indeed, the quantities used 
to construct this matrix can all be measured experimentally using the techniques developed 
by Ferris et al. |2l!]. What we have also seen here that was not visible in the analysis used 
by Strzys and Anglin is that the details of the approach to the final dynamical equilibrium 
state depend strongly on the initial quantum states. 



VII. RELAXATION TO EQUILIBRIUM 



The relaxation to equi 
as seen in, for example 



ibrium of closed quantum systems is an important topic of study. 
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23| . with a beautiful experiment by Kinoshita et al. [24] having 
shown that relaxation to equilibrium does not happen in a trapped one dimensional Bose gas. 
This was not unexpected for a one dimensional untrapped Bose gas with point interactions, 
which is known to be an integrable system, but it had been thought that practical features 
such as the harmonic trap and imperfectly point-like interactions would compromise the 
integrability and the system would relax. Our system is not integrable and is therefore 
able to equilibriate at zero temperature, without any interactions with a thermal cloud or 
other reservoir. However, to the best of our knowledge, there is as yet no consensus on the 
mechanism by which this happens. 
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For the thermahsation of quantum systems at finite temperature, one proposal is the 



eigenstate thermahsation hypothesis 
imphcitly contains a thermal state j 



TH), in which every eigenstate of the Hamiltonian 



25|. Srednicki, when introducing this hypothesis, 



26|, which is 



claimed that a necessary condition was the validity of Berry's conjecture 
expected to hold for systems which exhibit classical chaos in at least a large majority of 
the classical phase space. It is a simple matter to calculate effective Lyapunov exponents 
for the aj in this system, as long as we restrict the inital region of the phase space to that 
close to the initial conditions we have used for well occupations, and thus determine whether 
classical chaos could be present. A Lyapunov exponent for each well can be defined as 

h = li- -^T^. (25) 

where 

5a,{T) = \af\T)-af\T)l j = 1,2,3,4, (26) 

where is an initial condition slightly perturbed from In practice, we obviously 

cannot integrate the equations for infinite time, so we integrate the coupled GPE type 
equations over a reasonably long time and look at the development of 5aj{t) and hence 
Lj{t). What we found was that the system was stable for initial distributions close to equal 
numbers in each well, but that, for the initial numbers used in Fig. [121 it was unstable 
and therefore chaotic, meaning that it satisfies the criterion for Berry's conjecture to be 
applicable. The Lyapunov exponents as a function of time are shown in Fig. [H] for this 
unstable configuration. 

The fact that the classical equivalent of our system is chaotic for some initial conditions 
and stable for others suggests that it may provide a useful laboratory for the investigation 
of thermahsation in closed quantum systems. It also explains why we do not expect to see 
full collapses and revivals in the quantum system for arbitrary initial conditions, in contrast 
to those found for twin wells, where there are as many constants of the motion as there are 
equations of motion. Although it is prohibitively difficult to calculate either the eigenstates 
of the system or the full density matrix, we have developed practical alternatives which 
should be experimentally measurable. Further investigation of this model in terms of zero 
temperature thermahsation will be a subject of future study. 
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FIG. 14: (colour online) The Lyupanov exponent Li(t) for the same parameters as used in Fig. [THl 
The solid line has one atom added to well oi and one subtracted from 62 as the perturbation, while 
the dash-dotted line has a difference of two atoms. 

VIII. CONCLUSIONS 

We have used stochastic integration in the truncated Wigner representation to examine 
the four-mode Bose-Hubbard model proposed by Strzys and Anglin as a model for heat 
transfer, in parameter regimes not considered by the original authors. The inclusion of more 
quantum effects in our analysis shows that this system is not well described by linearisation 
of fluctuations about classical solutions, and moreover, that the initial quantum states used 
have a qualitative effect on the subsequent dynamics. The semi-classical analysis predicts 
that the system will exhibit First Law phenomenology, with continuing oscillations between 
different types of energy, and will not exhibit irreversible spontaneous processes which would 
result in an increase of the system entropy. Our quantum analysis suggests that the oscilla- 
tions will not necessarily be persistent and that an analogue of irreversible processes takes 
place which does lead to an increase in entropy, with the system therefore exhibiting both 
diffusive and oscillatory behaviours. 

We suggest that the model may not be a good one for heat transfer, for the reason that it 
is not at all obvious what the equivalent of temperature is for this system. As temperature 
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will flow from a hotter to colder object, it cannot be the number of atoms in a well, as the 
overall flow of atoms is at times from a less occupied to a more occupied well. However, the 
freedom of parameter choices in the initial conditions and the presence of both stable and 
unstable regions of the initial classical phase space suggest that it will be a highly tractable 
model for closed system quantum thermalisation. 
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